      SUBROUTINE E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,A,LA,IT,
     *                  RNMBST,RNM,IMPROV,ADIF,RES,PMAX,PINDEX)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988).
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE PRESID
C
C     CREATED 18 02 80.  UPDATED 14 05 80.  RELEASE 00/28
C
C     AUTHOR ... MAURICE G. COX.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEV.  FORMS PERFORMANCE INDICES AND RESIDUALS
C     FOR A POLYNOMIAL APPROXIMATION  P(X)  TO A SET OF DATA
C     WHICH MAY CONTAIN DERIVATIVE VALUES.  ALSO INDICATES
C     WHETHER THE POLYNOMIAL DEFINED BY THE SPECIFIED
C     COEFFICIENTS IS A BETTER APPROXIMATION THAN THE BEST
C     SO FAR.
C
C     INPUT PARAMETERS
C        WITHQ0   TRUE IF ZEROIZING POLYNOMIAL, ELSE FALSE
C        WITHPI   IF TRUE, PERFORMANCE INDICES PRODUCED UNDER
C                    ALL CIRCUMSTANCES.  OTHERWISE, PRODUCED
C                    ONLY IF  P(X)  IS AN IMPROVEMENT
C        M        NUMBER OF X-VALUES.  ALL DISTINCT
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C        X        X-VALUES.  NORMALIZED TO  (-1, 1)
C        N        NUMBER OF Y-VALUES
C        Y        VALUES AND DERIVATIVES OF DEPENDENT VARIABLE
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE
C        IMAX     ONE GREATER THAN LARGEST VALUE OF  IP
C        A        CHEBYSHEV COEFFICIENTS OF  P(X)
C        LA       DIMENSION OF  A  AND  ADIF.
C                    .GE. N + 1 IF WITHQ0 IS TRUE,
C                    .GE. N     OTHERWISE
C        IT       ITERATION NUMBER
C
C     INPUT/OUTPUT PARAMETERS
C        RNMBST   2-NORMS OF RESIDUALS CORRESPONDING TO THE
C                    BEST POLYNOMIAL SO FAR AND ITS DERIVATIVES
C
C     OUTPUT PARAMETERS
C        RNM      2-NORMS OF RESIDUALS CORRESPONDING TO
C                    P(X)  AND ITS DERIVATIVES
C        IMPROV   TRUE IF  P(X)  IS AN IMPROVEMENT, ELSE FALSE
C        ADIF     CHEBYSHEV COEFFICIENTS OF  (IMAX - 1)-ST
C                    DERIVATIVE OF  P(X)
C        RES      RESIDUALS CORRESPONDING TO Y-VALUES
C      * PMAX     LARGEST PERFORMANCE INDEX
C      * PINDEX   PERFORMANCE INDICES
C
C     NOTE.  THE PARAMETERS MARKED  *  ARE PROVIDED ONLY
C            IF  IMPROV  IS TRUE
C
C     KEY LOCAL VARIABLES
C        ASUMAX   FOR CURRENT VALUE OF  L,  MAXIMUM VALUE
C                    OVER  J = 0, 1, ..., L  OF SUM OF MODULI
C                    OF CHEBYSHEV COEFFICIENTS OF DERIVATIVE
C                    OF ORDER  J  OF  P(X)
C        EPS      RELATIVE MACHINE PRECISION
C        IY       LOCATION IN  Y  OF CURRENT SPECIFIED
C                    DERIVATIVE VALUE OF ORDER  L
C        NL       NUMBER OF SPECIFIED DERIVATIVE VALUES
C                    OF ORDER  L
C        NTERMS   NUMBER OF TERMS IN CHEBYSHEV SERIES
C                    REPRESENTATION OF CURRENT DERIVATIVE
C                    OF  P(X)
C        RES2NM   2-NORM OF THE RESIDUALS CORRESPONDING TO THE
C                    SPECIFIED DERIVATIVE VALUES OF ORDER  L - 1
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PMAX, XMAX, XMIN
      INTEGER           IMAX, IT, LA, M, N
      LOGICAL           IMPROV, WITHPI, WITHQ0
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LA), PINDEX(IMAX), RES(N),
     *                  RNM(IMAX), RNMBST(IMAX), X(M), Y(N)
      INTEGER           IP(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSRES, ASUM, ASUMAX, EPS, HALF, MLTPLR, ONE, P,
     *                  PLARGE, RESMAX, RSCALE, T, ZERO
      INTEGER           I, IA, IY, L, LM1, NL, NTERMS
      LOGICAL           IMP
C     .. External Functions ..
      DOUBLE PRECISION  X02AJF
      EXTERNAL          X02AJF
C     .. External Subroutines ..
      EXTERNAL          E02AHZ, E02AKZ
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SQRT
C     .. Data statements ..
      DATA              HALF, ZERO, ONE, MLTPLR/0.5D+0, 0.0D+0, 1.0D+0,
     *                  8.0D+0/
C     .. Executable Statements ..
      PMAX = ZERO
      EPS = X02AJF()
      NTERMS = N
      IF (WITHQ0) NTERMS = NTERMS + 1
      ASUMAX = ZERO
      DO 20 I = 1, NTERMS
         ADIF(I) = A(I)
   20 CONTINUE
      NTERMS = NTERMS + 1
      DO 120 L = 1, IMAX
         NTERMS = NTERMS - 1
C
C        DETERMINE SUM OF MODULI  ASUM  OF CHEBYSHEV COEFFICIENTS
C        OF DERIVATIVE OF ORDER  L - 1  OF  P(X),  AND UPDATE
C        ASUMAX
C
         ASUM = HALF*ABS(ADIF(1))
         IF (NTERMS.EQ.1) GO TO 60
         DO 40 IA = 2, NTERMS
            ASUM = ASUM + ABS(ADIF(IA))
   40    CONTINUE
   60    IF (ASUM.GT.ASUMAX) ASUMAX = ASUM
C
C        PINDEX(L)  IS USED TEMPORARILY TO HOLD THE
C        VALUE OF  ASUMAX
C
         PINDEX(L) = ASUMAX
         IY = L
         NL = 0
C
C        TO REDUCE THE POSSIBILITY OF UNDERFLOW AND OVERFLOW,
C        COMPUTE  RES2NM  AS  RESMAX*SQRT(RSCALE),  WHERE
C        RESMAX  AND  RSCALE  ARE UPDATED FOR EACH RESIDUAL
C        CORRESPONDING TO A SPECIFIED DERIVATIVE VALUE OF
C        ORDER  L - 1.  AT ANY STAGE,  RESMAX  HOLDS THE
C        MODULUS OF THE RESIDUAL OF MAXIMUM MAGNITUDE SO FAR,
C        AND  RSCALE  THE SUM SO FAR OF THE SQUARES OF THE
C        RESIDUALS, EACH SCALED BY  RESMAX.
C
         RESMAX = ZERO
         RSCALE = ONE
         DO 100 I = 1, M
C
C           SKIP IF NO DERIVATIVE VALUE OF ORDER  L - 1  IS
C           SPECIFIED AT  I-TH  X-VALUE
C
            IF (IP(I)+1.LT.L) GO TO 80
            NL = NL + 1
C
C           EVALUATE  P,  THE  (L - 1)-ST  DERIVATIVE OF  P(X)
C           AT  X = X(I)
C
            CALL E02AKZ(NTERMS,ADIF,1,LA,X(I),P)
C
C           SAVE RESIDUAL CORRESPONDING TO THIS VALUE
C
            IF (WITHQ0) RES(IY) = -P
            IF ( .NOT. WITHQ0) RES(IY) = Y(IY) - P
C
C           UPDATE  RESMAX  AND  RSCALE
C
            ABSRES = ABS(RES(IY))
            IF (ABSRES.NE.ZERO .AND. ABSRES.LE.RESMAX) RSCALE = RSCALE +
     *          (ABSRES/RESMAX)**2
            IF (ABSRES.NE.ZERO .AND. ABSRES.GT.RESMAX)
     *          RSCALE = RSCALE*(RESMAX/ABSRES)**2 + ONE
            IF (ABSRES.NE.ZERO .AND. ABSRES.GT.RESMAX) RESMAX = ABSRES
   80       IY = IY + IP(I) + 1
  100    CONTINUE
         RNM(L) = RESMAX*SQRT(RSCALE)
C
C        FORM COEFFICIENTS IN CHEBYSHEV SERIES REPRESENTATION
C        OF  L-TH  DERIVATIVE OF  P(X)  FROM THOSE OF
C        (L - 1)-ST DERIVATIVE
C
         IF (L.LT.IMAX) CALL E02AHZ(NTERMS,XMIN,XMAX,ADIF,1,LA,T,ADIF,1,
     *                              LA)
  120 CONTINUE
C
C     IF NOT ON ZERO-TH ITERATION,
C     DETECT WHETHER THERE HAS BEEN AN IMPROVEMENT ...
C
      IMP = (IT.EQ.0)
      IF (IMP) GO TO 160
      DO 140 L = 1, IMAX
         IF (RNM(L).LT.RNMBST(L)) IMP = .TRUE.
  140 CONTINUE
      IF ( .NOT. (IMP .OR. WITHPI)) GO TO 220
C
C     ... AND, IF SO, OR IF ZERO-TH ITERATION, OR IF
C     THEY ARE REQUIRED ANYWAY, FORM THE PERFORMANCE
C     INDICES CORRESPONDING TO THE IMPROVED
C     APPROXIMATION
C
  160 PLARGE = ZERO
      DO 200 L = 1, IMAX
C
C        NL  IS THE NUMBER OF DERIVATIVE VALUES OF ORDER  L - 1
C
         NL = 0
         LM1 = L - 1
         DO 180 I = 1, M
            IF (IP(I).GE.LM1) NL = NL + 1
  180    CONTINUE
         T = NL
C
C        RETRIEVE THE VALUE OF ASUMAX
C
         ASUMAX = PINDEX(L)
         IF (ASUMAX.NE.ZERO) PINDEX(L) = RNM(L)
     *       /(MLTPLR*EPS*ASUMAX*SQRT(T))
         IF (ASUMAX.EQ.ZERO) PINDEX(L) = ZERO
         IF (PINDEX(L).GT.PLARGE) PLARGE = PINDEX(L)
  200 CONTINUE
      PMAX = PLARGE
  220 IMPROV = IMP
      RETURN
C
C     END E01AEV
C
      END
